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Abstract 

We present the first fifth-order, semi-discrete central-upwind method for ap- 
proximating solutions of multi-dimensional Hamilton-Jacobi equations. Unlike 
most of the commonly used high-order upwind schemes, our scheme is formulated 
as a Godunov-type scheme. The scheme is based on the fluxes of kurganov- 
Tadmor and Kurganov-Tadmor-Petrova, and is derived for an arbitrary number 
of space dimensions. A theorem establishing the monotonicity of these fluxes 
is provided. The spacial discretization is based on a weighted essentially non- 
oscillatory reconstruction of the derivative. The accuracy and stability properties 
of our scheme are demonstrated in a variety of examples. A comparison between 
our method and other fifth-order schemes for Hamilton-Jacobi equations shows 
that our method exhibits smaller errors without any increase in the complexity of 
the computations. 


Key words. Hamilton-Jacobi equations, central schemes, semi-discrete schemes, high 
order, WENO, CWENO, monotone fluxes. 
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1 Introduction 


We are interested in approximating solutions of multi-dimensional Hamilton-Jacobi (HJ) 
equations of the form 

4> t + = 0, x = (x!,...x d ) € ( L1 ) 


where 6 = 4>{x,t), and the Hamiltonian, H. depends on Vd> and possibly on x and U 
Solutions of (1.1) develop discontinuous derivatives even for smooth initial data, lhis 
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loss of regularity presents difficulties both in the analysis of these equations as well as 
in numerically approximating their solutions. Significant advances in the theoretical 
understanding of the HJ equations were achieved in the last two decades. Most notable 
is the introduction of the so-called “viscosity solution” which provides a consistent def- 
inition of a weak solution of (1.1) past the time where the discontinuities develop. See 
[2, 7, 8, 9, 10, 14, 25, 26] and the references therein. 

In spite of the large number of applications for H J equations, there has been very 
little activity in numerically approximating their solutions. This is surprising in partic- 
ular given the equivalence between the HJ equations and hyperbolic conservation laws, 
and the flourishing field of numerical methods for conservation laws. Converging first 
order methods for the HJ equations were introduced by Souganidis in [33]. High order 
upwind methods were introduced by Osher, Sethian and Shu in [31, 32]. The schemes m 
[31, 32] were based on an “essentially non-oscillatory” (ENO) reconstruction by Harten 
[13] and a monotone numerical flux. A more compact upwind scheme which is based on 
a weighted ENO (WENO) reconstruction is due to Jiang and Peng [15]. WENO recon- 
structions were originally introduced in the context of numerical schemes for hyperbolic 
conservation laws in [16, 29]. They increase the order of accuracy by using wider stencils 
in smooth regions while automatically switching into one-sided stencils in regions t at 
include singularities. All these reconstructions include nonlinear limiters in order to 
control the spurious oscillations that might develop in the solution. For extensions to 

unstructured grids see [1, 34]. , . 4 , . 

A class of Godunov-type approximations for HJ equations was recently introduced 
by Lin and Tadmor in [27, 28]. Their first- and second-order central schemes were 
based on the first-order Lax- Friedrichs scheme [11] and the second-order Nessyahu- 
Tadmor scheme [30] for approximating solutions of hyperbolic conservation laws. Cen- 
tral schemes incorporate internal averaging over discontinuities and hence they do not 
require Riemann solvers. Moreover, systems can be solved without a characteristic de- 
composition and this makes central schemes simple, robust, and particularly suitable 
for treating complex geometries. We developed in [4] an efficient version of the central 
schemes of [27, 28] for multi-dimensional HJ equations. Our first- and second-order, 
non-oscillatorv, non-staggered schemes were designed to scale well with an increasing 
dimension. Efficiency was obtained by carefully choosing the location of the evolution 
points and by using a one-dimensional projection step. In [5, 6] we introduced third- 
and fifth-order fully-discrete central schemes, which were the first central schemes for 
HJ equations of accuracy greater than two. High-order accuracy was obtained using 
a suitable high-order WENO-type reconstruction. We would like to note that ENO 
and WENO interpolants were already used in central schemes for conservation laws 

[3, 22, 23, 24]. 

One way to improve the above schemes [4, 5, 6, 27, 28] is to use semi-discrete meth 
ods to reduce the numerical dissipation. In principle, one expects to obtain a semi- 
discrete scheme from a fully-discrete scheme in the limit as At -*• 0. Unfortunately, 
the fully-discrete schemes in [4, 5, 6, 27, 28] blow up in that limit. A different strategy 
is to consider at every grid point more precise information regarding the local speed 
of propagation, which can then be used to develop a different class of fully-discrete 
approximations that do enjoy a semi-discrete limit. An estimate of the local speed of 
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propagation at every grid point can then be used to determine new points where the 
solution is evolved to the next time step. The distance of these evolution points from the 
original grid points is proportional to the time step At, and hence it is possible to ob- 
tain a semi-discrete scheme in the limit At -> 0. This strategy was first used to develop 
semi-discrete central schemes for hyperbolic conservation laws: a second-order metho 
was developed by Kurganov and Tadmor in [21], and a third-order method by Kurganov 
and Levy in [18]. Semi-discrete schemes for HJ equations were then introduced m [ J, 
and further improved in [19] by utilizing a more accurate estimate of the local speed of 
propagation, hence reducing numerical dissipation. We would like to stress that bot 

schemes [19, 20] are only second-ox der accurate. 

In this paper we present fifth-order, semi-discrete, Godunov-type, central schemes 
for HJ equations. These are the first high-order semi-discrete central schemes for HJ 
equations 1 These schemes are generated by considering a general formulation of semi- 
discrete schemes along the lines of [19, 20], and augmenting it with an appropriate 
high-order WENO-type reconstruction. 

The structure of this paper is as follows. In Section 2 we develop a one-dimensiona 
fifth-order semi-discrete scheme. In the semi-discrete limit, A t — * 0, the t or er 
WENO interpolant we obtain turns to be identical to the one used m upwind methods 
[15] . The flux we use is the Kurganov-Noelle-Petrova flux [19], or a variant o t e 
simpler Kurganov-Tadmor flux [20]. We state a theorem establishing the monotonicity 
of these fluxes, the proof of which is left to an appendix. We observe that for the 
one-dimensional linear advection, our method boils down to an upwind scheme with 
a Lax-Friedrichs flux. In Section 3 we generalize the method to an arbitrary number 
of space dimensions, writing out the schemes explicitly for two and three dimensions 
in Section 3.2. We conclude in Section 4 with several numerical examples m one, two 
and three space dimensions that confirm the expected order of accuracy as well as the 
high-resolution nature of our scheme. We compare the results of these numerical tests 
with our fully-discrete fifth-order scheme [6] and with the scheme of Jiang and Peng [ T 
Our numerical results show that the new method we present in this paper has stability 
properties that are equivalent to those of [15]. The relative V errors we obtain m all 
our simulations are consistently smaller than those in [15], in some cases as much as an 
order of magnitude smaller. 

Acknowledgment: The work of D. Levy was supported in part by the National Science 
Foundation under Career Grant No. DMS-0133511. 


2 A One-Dimensional Scheme 

2.1 Semi-Discrete Central Schemes for HJ Equations 

Consider the one-dimensional HJ equation of the form 

Mx,t) + H (<j> x ) = 0, x€R. 


1 high-order is assumed to be an order greater than two. 
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Figure 2.1: The regions considered in Godunov-type central schemes. The ± solution is 
evolved at x, ± af At. The solution <pf +x is obtained by averaging <p(x t ± a , A t,t n+ ). 


We are interested in approximating solutions of (2.1) subject to the initial data 
0) = (f>o{x). We briefly review the construction of semi-discrete central schemes for 
(2.1) presented in [19] (see also [20]). For simplicity we assume a uniform grid grid in 
space and time with mesh spacings. Ax and At, respectively. Denote the grid points 
ky x . __ t n — nAt. Let ip™ denote the approximate value of and at a 

fixed time t n let <p\ denote the approximate value of the spatial derivative (p x (x l ,t ). We 
define the forward and backward differences A + (p t := <p i+ 1 - pi and A pi := p x - Pi-i- 
Assume that the approximate solution at time t n , pf is given, and that a continuous 
piecewise-polynomial interpolant p(x, t n ) w r as reconstructed from yh ■ The construction 
of <p(x t n ) will be addressed below. At every grid point x r we then estimate the maximal 
speed of propagation to left, of, and to the right, af . For a convex Hamiltonian, these 
one-sided local speeds of propagation are estimated bj 


af = max {H' (^' ) , H' (<^( + ) , 0} , a t 


min {H' [ip\ ) , H' (v^[ + ) . 0} | • (2-2) 


Here, ip’± are the one-sided derivatives, defined as 
:= lim <p x (x, ± af At, t n ) . 

At—0 

Remark. Our sign convention in the definition of af in (2.2) differs from [19]. This 
choice of signs simplifies the derivation of the scheme in the multi-dimensional setup. 

We evolve (p according to (2.1) at the evolution points x, ± a. At. The time step, 
At, is chosen so that the reconstruction is smooth at these points (see Figure 2.1). A 
Taylor expansion in time of dp (x, ± af At, f n+1 ) results with 

dp (x, ± afAt, t n+l ) =<p(x,± af At, t n ) - AtH (< p x (x, ± of At, t n ))+0 (At 2 ) .(2.3) 

A weighted average is then used to re-project dp (x, ± afAt, t n+v ) onto the original grid 
point x i) 


_n+l 

S °i - 


af + a, 


-dp (x, - a i At, t n+1 ) + 


af + a, 


'W 


5 (x t + af At. t n+1 ) . 


(2.4) 
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A fully discrete central scheme is finally obtained by substituting (2.3) into (2.4) 

n+i __ (<^( Xi -a-Af,t n ) - AtH (p x {xi-a~At,t n ))) (2.5) 

af + a* 

-I — — — (p (x^ + af At, t n ) — A tH {fp x (x^ + a, At, t ))) ■ 

af + a { 

Utilizing the Taylor expansion p ( x< ± af At,t n ) = p(xi,t n ) ±afAtp , and assuming 
that p satisfies the interpolation requirements pf = p(x,,t n ), equation (..5) can 
rewritten as 


Vi 


,n+l _ , n n 


pf + At [vt ~ Vx] 


(2.6) 


d: 4*" Cl 


At 


af + a, 


[a { H (p x (xi + afAt.,t n )) + a fH (Vx a i ))] ' 


Here <5± denotes the one-sided reconstruction of the derivative at x t . A Godunov-type 
semi-discrete method for approximating the solution of (2.1) is obtained by taking the 
limit At — v 0 in (2.6) (see [19, Eq. (3.44)]) 


— pi (t) = --qr- — r [a, H {v't) + a ? H {v\ )] + 
dt 


af + O; 


(p[ + - V?) • M 

a f + a i 


Even though the flux on the right hand side of 2.7 was originally presented in [ ], 

Kurganov et. al. did not investigate its properties. We now state a theorem establishing 
the monotonicity of this flux. The proof is left to the appendix. 


Theorem 2.1 Assume that H € is convex. Then 

1 


H knp [u + ,u~) 


a 


■H ( u + )+a + H ( u )] 


a + a 


( u + - u ) , 


a + + a- L ~~' 7 ' a+ + a- 

is a monotone flux, that is H KNP is a non-increasing function of u + and a non- 
decreasing function of u . 


Remarks. 

1 The derivation of the semi-discrete scheme (2.7) does not depend on choice of 
interpolants p, so long as they are smooth at x, ± af Af during the time interval 
[t n ,t n+1 l. The spatial order of accuracy of the scheme is determined by the 
accuracy of the reconstruction of p as well as on the accuracy of the ODE solver 
used to solve (2.7). A suitable high-order reconstruction will be presented m 
Section 2.2 below. To be precise, the scheme (2.7) does not require the 
construction of the interpolant p. All that is required is a reconstruction of the 
point- values of the derivatives p*. 
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Figure 2.2: The three interpolants used for the fifth-order reconstruction + t + - In this 
example, because of the large gradient between x i+i and x,+ 2 , the interpolant <?- Wl1 
have the strongest contribution to the CWENO reconstruction at x l+T . 


2. In order to economize on storage space, and sometimes also reduce the 
computations, it is possible to replace a t + and a- with 
a t = max {| H' (^)h \ H> (<^D | } • In this case ’ ( 2J ) becomeS 

J t <Pi (t) = ~\ [ H (<^ + ) + H faD] + J ( 

This simpler formulation was presented by Kurganov and Tadmor in [20, Eq. 
(4.10)]. We denote the right hand side of (2.8) by -H KT )- 

As an immediate consequence of Theorem 2.1, we have 

Corollary 2.1 If H € C 2 is convex, then H KT {u + ,u~) is a monotone flux. 


2.2 A Fifth-Order Scheme 

In order to obtain a fifth-order scheme from the general semi-discrete formulation (2.7) 
we need a fifth-order approximation of the derivative *>' and a suitable ODE solver. A 
central-upwind interpolant at Xi starts with a central interpolant defined either on t e 
interval [xi,x i+ 1 ] for a right-biased reconstruction, or [xi-i,x<] for a left- biased recon- 
struction. This central interpolant is then evaluated at the location x t+T x, + r Ax, 
where r is a parameter introduced for notational convenience. For the semi-discrete 
scheme (2.7) we take r = afAt for the right-biased interpolant, and r = o. At for the 

left-biased interpolant. , . . A , , 

For the right-biased interpolant at x i+T (r = afAt), we use three cubic interpolants 

,J~ k = l, 2, 3, defined on the stencil {xi-3+fc; . . . , x l+ k )• Here 

r k.l + T ? 3 


^i + ,+r 


^2 + i+r 


1 


6Ax 


(1 - 3r 2 )<p,_2 + 3(— 2 + 2 r A 3r 2 )^. 


i— 1 


+3(1 - 4r - 3 r 2 )ipi + (2 + 6r + 3r 2 )<^+i | 

— - — (— 2 + 6t — 3T 2 )<f?j-_i + 3(— 1 — 4r + 3r )(p t 
6Ax L 

+3(2 + 2r - 3r 2 )^+i + (-1 + 3r 2 )v> i+2 


(2.9) 
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VZ + T = 


6Ax 


(-11 + 12t - 3r 2 )(y?i + 3(6 - 10r + 3 r 2 )<p I+1 


+3(-3 + 8r - 3r 2 )ifi+2 + (2 - 6r + 3 t 2 )^ +3 J 


A straightforward computation shows that V/c, <P fc + , +T — a x V ( x i+t) + 0 ((Ax) ). Also 
the following linear combination is a fifth-order approximation of ip t + 


Vitr 


— ^2 C k<Pk,i+T ip + O ((Ax) 5 ) , 


fc= 1 


provided that the constants Ck are taken as 

! 15t 2 + IOt - 6 - 120r 3 + 120r 4 


Cl 20 


c? = — — 


3t 2 — 1 

1 720r 6 - 1080r 5 + 660r 4 + 6 0 t 3 - 81r 2 - 64r + 24 
20 (3 t 2 - 1) (2 - 6r + 3r 2 ) 

1 — 15r 2 + 4 + 120r 4 


° 3 20 2 — 6r + 3r 2 

In the limit r->0, y'Z := lim T _ o^ +T — al 
1 


, 0 '+ 


, r f + 

V2,i 


6Ax 


§f (x,) + 0((Ax) 3 ), with 

(ipi-2 — 6f^i-i + 3 <fli + 2(pi + i) , 


— ( — 2 / >Xi — X — 3<p t 6^j-|-X ^*+ 2 ) 

6Ax 


m'+ = — (-11^, + 18^+1 - 9^i+2 + 2 { Pi+Z ' 

6 Ax 

A right-biased fifth-order interpolant at x t is therefore given by 


r = = I:* w + 0 «■ * x)6) • 


Vi 


10 


10 


(2.10) 


By symmetry, for the left-biased interpolant (r = a" At) we use three cubic inter- 
polants V fc "+T> * = 1,2, 3, this time defined on the stencil {x t _ 4 +fc, • • • . x i- i+* /■ In 1 e 
limit r - 0, Vm == ^ = £(*<) + ^((Ax) 3 ), where 


Vi.i = 


V2.» — 


= 


1 


6Ax 

1 

6Ax 

1 

6Ax 


(2<Pi_3 — 2 + 18(p,-l — HV») I 

(— + 6i/?j_x — 3<Pi 2(pi+i) . 

(2<Pi_l + 3ip, - 6iy2i + l + <^*+2) • 


In this case 


V,' 




a 


10 ' 


10 


( 2 . 11 ) 
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In order to suppress spurious oscillations, the coefficients in S are replaced by nonlinear 
weights, which are set as to preserve the order of accuracy of the reconstruction in 
smooth regions while automatically switching to the appropriate stencil in regions that 
contain discontinuities. To this end we define the convex combination 


<Pi 


-E 


-4- j~ 


h = 1 


E 

fc=i 


= 1. 


( 2 . 12 ) 


In smooth regions ~ « cf — c 3 — 3/10, ~ w l { ~ c 3 c i 1/10 and 

± ~ ^ = 3/5 so the order is of the order O ((Ax) 5 ). When the stencil supporting 
contains a discontinuity, the weight of the more oscillatory polynomial should vanish. 
Following [16, 29], these requirements are met by setting 


UJ 


= 


fc,i 


± > 




CJ 


Ar.i 




(2.13) 


where it l € {1,2.3}. We choose e as 1(T 6 to prevent the denominator in (2.13) from 
vanishing, and set p = 2 (see [16]). The smoothness measures should be large when 
V is nearly singular. Following [16], we take S± to be the sum of the L 2 -norms of the 
derivatives on the stencil supporting tp' k i . If we approximate the first derivative at x, 
by A + ^/Ai, the second derivative by A + A _ ^ I±r /(Ax) 2 , and define the smoothness 


measure 


S^s] = Ax^ (^A + y>i+j) + A;r ^ (ax 2A+A ^ +i ) ’ 


(2.14) 


then for the right-biased interpolant we have Sf t =^5 t [-2, 0], S£ t - S, [ 1,1] and 
5+ = Si [0, 2]. For left-biased interpolant we have = 5, [-3, -1], S 2<i = S t [-2, 0] 
and = 5, [-1, 1]. We use the notation 

ip r± = reconstruct.^' (±, S ) , (2-15) 


to denote the computation of the array {y?'*} for all i from data y? n at time t n , as given 

by (2-12). , . , , 

The following algorithm summarizes our fifth-order semi-discrete algorithm tor ap- 
proximatinng the solution of (2.7). The time integration is performed with a fourth-order 
strong stability preserving (SSP) Runge-Kutta scheme [12]. 

Algorithm 2.1 Let F (<f'~ , y?' + ) denote the right hand side of (2.7). Then at each grid 
node i, 


tp'tf = reconstruct S (- , y>”) , S = reconstruct S (+, S) 


^) = y, n + i AtF {&,<&) 

if'- = reconstructs (-, yj (1) ) , v’i = reconstructs (+, y? (1) ) 
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(2 ) = _649_ 

^ 1600 v 


10890423 aif (*,#) + + H^f (^n 


F2 


25193600 

n(2) 


/- = reconstruct jpf (-, F (2) ) , F2 + = reconstruct jpf (+, F (2) ) 


,(3) 




F3 


53989 


2500000 
5121 


20000 


^ _ 5000000 AtF ,V?0 ^ + 20000000^ 

10000 v 2 ' 


, ,, V Ui3 

A^F (Fi >V?i ) + 32000^ 


23619 (2) + 


: reconstruct jpf (-, <p (3) ) , & = reconstruct. jft (+, F (3) ) 


<P 


n + 1 


V + ^AtF (<Fo',Fo + ) + 
o 1U 

i(^ (3) + jj;AtF (9?3 _ ,^3 + ) 


6127 

30000 


F (1) + UtF (y/f,^) + 

0 


_7873_ (2) 
30000 


Remarks. 

1. The smoothness measures (2.14) are not the same as those used in [15, 16). 
There a different normalization of the derivatives was used. Our smoothness 
measures are approximations to the sum of the L 2 - norms of the first and second 
derivatives of the interpolant on a stencil. In the cases we tested our smoothness 
measures produced comparable or smaller errors when compared with [15]. We 
include a comparison between the results obtained with both forms of the 
smoothness measures in Section 4.1.1. 

2. From obvious reasons, the interpolant (2.10) is identical to the one used in the 
upwind method of [15]. As far as the scheme itself is concerned, there is some 
degree of similarity between the semi-discrete central scheme and upwind 
schemes. It is important to note that for linear advection problems they boil 
down to the same scheme. Indeed, if H(s) = s , then H' — 1. Hence a- — 1, 
a~ = 0 and equation (2.7) becomes 


d_ 

dt 


Fj (*) = 




a j 


(2.16) 


Solving (2.16) is equivalent to solving 


!*,,(() = -tf LF W*, ¥>;-), 

with the Lax-Friedrichs flux 

H LF (f[ + , ip ’~) = H Q (Fi + + Fi _ )) - g ^ ^ = ifl ' 
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For this reason the schemes in [19] are called “central-upwind schemes”. Even in 
this case of a linear advection problem, there still are some differences between 
our scheme (Algorithm 2.1) and the scheme in [15]: the ODE solvers and the 
smoothness measures are different. For more complicated Hamiltonians the 
semi-discrete scheme (2.7) is different than the scheme in [15]. A comparison 
between numerical results obtained with both schemes can be found in Section 4. 

3 One can easily create a third-order semi-discrete central scheme from the general 
one-dimensional formulation (2.7) by using a less accurate ODE solver and a 
third-order interpolant. Indeed, a third-order version of the right-biased 
(derivatives of the) interpolants can be written as a combination of two 
polynomials, ip' u with j = 1,2, that are constructed on the stencil 
[ X]+i _ 2 , . . . , x J+i } (compare with (2.9)). A straightforward computation shows 

that 




V 3 2,i 


lim — — 
t— » o Ax 


1 


- + T ) Ipi - 1 + ( — 2d) p>i + [ 0 -r T J ‘fi+l 


2Ax 


{‘Pi+l *Pi- 1) i 


lim — — 
t — > o Ax 


-- + T ] Ifii + (2 - 2r ) + I + T ) ^>+2 


—7— (-3 ifii + 4<^j+i — <fi+ 2) 
2Ax 


satisfies <*+.) /dx + O ((A for t * 0 and j = 1, 2. The combination 

= >!3 


1 2 - 6r + 3r 2 , 1 — 1 + 3r , 

+ 3 _i + 2t 


-1 + 2r 


2 , 1 , 
= 3^ + 3^ 2 ’ 1 


satisfies ^ = dp (x i+ r) /dx + O ((Ax) 3 ) . The left-biased interpolants can be 
easily derived by symmetry considerations. 


3 Multi-Dimensional Schemes 


3.1 A General Multi-Dimensional Scheme 

Consider the d-dimensional HJ equation of the form 
d> t + H{V4>) = 0, x= (x (1) ,...,x (<i) ) € R d , 


(3.1) 


subject to the initial data <p(x,t — 0) — <po(x). 

For simplicity we assume a uniform grid in space Ax - = ■■■ = Ax - = Ax. We =>e 
a = (ai,a z ,...,a d )e Z d , and let x Q = Ax<5, such that the fc-th coordinate of x a equals 
X W _ Axa (fe \ VI < k < d. For example, in the conventional three-dimensional notation 
with indices i, j and k and components (x,y,z), a = (t,j,k) and x Q = (xi,yj,z k ). 
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Similarly to the one-dimensional setup, <p” will denote the approximation of <p(x a ,t ), 
and for a fixed time V<p Q will denote the approximation of V<£ at x Q . 

Given x a , we define the volume = (g>ti ~ + if] ’ and estimate the 

local speeds of propagation aj. For example, for a convex Hamiltonian these speeds in 
the coordinate direction k are given by 


a ( fe )+ = max 
° c a 


dH_ 

dxW 


(V<p a ),0 



(3.2) 


Let p = (p (1 \ . . . ,p (d) ) denote the multi-index with components p {k) € {+, -}, Vfc. We 
also denote the index opposite to p by p, i.e. p = -p = (— P (1 \ • ■ • > — P (d) )> assuming the 
standard algebraic operations between elements in Z. For any given p we define a vector 
that encodes the maximum estimated speed of propagation in all coordinate directions 

at x a as 


Z = (p (1) ai! 




■ 


(3-3) 


We then denote by x Q+p the position x Q + v p a A t, and denote the approximation of 4> at 
%a+p by <Pa+p- 


For example, if of = 3 and p — (+, — , +), then v p — , 

(-aL 1)_ ,aL 2)+ , -aL 3)_ ) ■ In this case 

< p q+ „ = v (& + <£ )+ a *, 4 2) - < 4 2)_ + a « )+ A *) ' 


(1)+ -aL 2) -,aL 3)+ ) 


and vf! = 


Similarly to the one-dimensional case, we assume that the approximate solution at 
time t n , <p” is given, and that a continuous piecewise-polynomial interpolant <p(x, t n ) 
was reconstructed from <p2- The construction of cp(x 7 t n ) will be addressed below. The 
interpolant <p(x, t n ) is then evolved to the next time step, t n+1 , at the points x a +p, which 
are located away from the propagating discontinuities assuming that the time step At 
is sufficiently small. According to (3.1), to first order in time, this evolution is given by 
the Taylor expansion 

<p (f a+p> t" +1 ) = <p(x a+p ,t n ) - A iff (Vp(x a+p ,n) + O {At 2 ) , 


where Vp {x a+p , t n ) is an approximation of the derivative V0 at x a+p . 

A fully discrete central scheme can then be constructed by computing a weighted 
average of the evolved solution p{x Q+p ,t n+l ) for all values of p (compare with the one-^ 
dimensional case (2.4)). The volume of the d-cube enclosed by x a + p for all values of p 
divided by At is 


v; = n W )+ + a « )_ ) • 

fc-i 

For a given p, the volume enclosed by the corners x a+p and x a divided by At is given 
by the product of the components of v p 

fc=l 
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Figure 3.1: A two-dimensional example of the objects associated with the location x a + p 
for p = (+, The complement location x a+p is shown, as well as the volumes |F£| and 
|u£|. The thick rectangle encloses the volume V a . 


Clearly, 1^1 = K>- See Figure 3.1 for a sketch of the two-dimensional setup. An 
approximation of the solution <p" +1 * s then obtained by averaging over all (p (x a+p , t + ). 
Each term corresponding to a particular p is weighted by the diagonally opposite volume 
r f | , divided by V a . Hence 


< +1 = 

Q P 

= ow n - AtH (V* ( x Q+P , r))i . 

P 


(3.4) 


We now use a Taylor expansion in space 

<p (x Q+p , t n ) = (p (x a ,t n ) + A tt£ • (x a+p , n + O (At 2 ) , 

where Vp(x Q+p ,t n ) is the evaluation of the gradient at x a associated with the recon- 
struction at x n+p in the appropriate volume. Hence (3.4) can be written as 


vr 1 


= vi + ttYI i^i k- v^(x a+P A n ) -H(vcp(x Q+p .,n)\ . 


V, 


In the limit At -> 0 we obtain our first form of the semi-discrete d-dimensional scheme: 

E i«S(*)i k • ^ a+P (t) - h (vp a+P m ( 3 - 5 ) 


± m__L 

dt^ j V a (t) 
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[v p a -V(p a+p {t)-H{Vy a+p {t))]. 


To obtain a simpler formula, we let, for p (k) = ±, = lim A t-.o d‘p>{x a+p )ldx , the 

k-th component of V<£ (x a+p ). Such a limit makes sense assuming that the reconstruction 
of the derivatives is done direction-by-direction. Then the first sum on the RHS of (3.5) 


becomes 


r E E 11 (*»+,) 

* a 


k= 1 p j = 1 

d 


t e e n 

* Ol 


k=l p j - 1 
d 


£ n «■ 


0)p (j) 


fc=l 


p 


H a 

a — i f * 

E tta tt-or / 4- — \ '^P — 

(fcV+ (/c)- 


/c=l 


ttL fc)+ + 


n U*- ( 


+ <£’“) 


11 a<f )+ a<.‘»- 


E cta 

+ a ^ 

tt Q T U,a 


r (^iffc) ^i(d) • 


This gives the semi-discrete d-dimensional scheme 
(f\ = - 

Ha 


d %<<) = -4 E Ki H ( v ^ +f ) + E 


d ai k)+ ai k) - 


dt 


a? )+ + og>- 


{'•pt(h) - V x w) ■ (3- 3 * * 6 ) 


iieraarfcs. 

1. The d-dimensional semi-discrete scheme (3.5) is valid for any reconstruction of 
V; p, including reconstructions defined on d-dimensional stencils (for 
two-dimensional examples see [6]). In contrast, (3.6) is valid only for 
dimension- by-dimension reconstructions such as those described in Section 3.3 
below. These dimension-by-dimension reconstructions are natural in the 
semi-discrete setting, as they significantly simplify the form of the scheme. 

2. As in the one-dimensional case, (3.5) and (3.6) are independent of the order of 
the reconstruction. First- and second order- reconstructions can be found, e.g ., 
in [19]. In Section 3.3 we develop a fifth-order dimension-bv-dimension 
reconstruction following the one-dimensional reconstruction of Section 2.2. 

3. Proof of the monotonicity of the flux approximation in (3.6) can be obtained via 

the method of proof of theorem 2.1 applied to each component. This becomes 

particularly transparent when (3.6) is written out as in Section 3.2 below. Such a 

proof cannot directly use the definitions aL fe)+ = maxc Q CV<p a ) , 0} etc., 


14 


S. Bryson and D. Levy 


where the maximum is taken over the spatial domain C a (see (3.2)). We must 
translate this definition into a maximum over the range of function values. For 
example, in two dimensions, we define 


a + = max {H x (u, v) ,0} , a 

„€/( u~,u+) 

C<v<V 


min {H x ( u , v ) , 0} 

( u + ) 

C<v<T> 


= max 

A< u<B 


{Hy(u,v),0}, b- = 


min {H y (u,v), 0} 

A<u<B 


where [. A , B] is the range of u and [C, V] is the range of v. With such a choice of 
a and 6 (and similarly in more than two space dimensions) the multi-dimensional 
flux approximation is monotone. 

4. The limit lim A t-o V<p(x a+P ) depends on p. Its value is estimated from the 
reconstruction that corresponds to p. 

5. When oi fc)+ and )_ are replaced by a ( n k) = max Co {\$Tj (Vy3 0 )|}, 

y a _ 2 d Y[ d k , a (k) = 2 d \v?\. In this case the semi-discrete scheme (3.6) becomes 

4v.(f) - -4 E H + \ E «« - aw) ■ (37) 

Z O k= 1 


A simpler one- and two-dimensional version of (3.7) was presented in [20] with a 
less accurate estimate of the local speed of propagation, a = max* a . 

6. In practice, the speeds of propagation is estimated from the reconstruction of 


a (fc)+ _ max 
P 


dH 

dx ^ k ) 


(V 7 j 0 



min 

p 


dH 

dx( k ) 


( V (jpa+p 



3.2 Two and Three Dimensional Schemes 

For convenience, we write out (3.6) in two and three dimensions. In two dimensions, we 
let a = (i,j) with coordinate notation (x x .y } ), and let the local speeds of propagation 

be (afj, bfj) := aL 2)± ) . Explicitly, for convex H we use the estimates 

a i"j = max { H x ,<Py) j 0 } 

Kj = max {H y (<p x ,ipy) ,0} 


a iJ = 


Kj = 


min {H x ,0} 

min {H y (<p * , <f% ) , 0} 


(3.8) 
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where the max and min are taken over all permutations of ±. Then (3.6) becomes 
(suppressing the indices i, j) 


dp 

dt 


a + a / + _\ , 

: (a + + a~) ^ ^ (b+ + b~) 

a~b~H (<p+, <p+) + a + b~H (pZ, pj) + a~5 + # (pi, Pi) + Q+6+ # 

(a+ + ar) (6+ + 6 _ ) 

If we replace aT and a~j by a t j = max± | { H x (pi , ) } | and similarly b iy and 

by b i: = max± \{H y (pi,pf)}\, then (3.9) can be further simplified to (compare with 
[20, Eq. (5.10)]) 


b + b~ 


-ZT « - Py ) 


(3.9) 


d J^ = a J±{vt-vZ) + h f{v + y -<p- y ) 

“ [H (<pt, <Py) + H + H ((ft , ’P~y ) + H (<PZ , Py )] ■ 


(3.10) 


In three dimensions, we let ol = with coordinate notation (xi, yj 7 Zfc), and let 

the local speeds of propagation be (a± jk , bf Jk , cf jk ) := (a2 ]± , aij , al ) j. Thus a and 
£> are the obvious generalization of (3.8), and c is estimated as 

ct jk = max {H : (pi, <py,<pf) > °} > c 7j,k = m ^ n i Hz fax ’ * ^ 2 ) > °) 


Then the semi-discrete scheme becomes (suppressing the indices i.j, k) 

1 


dp 

dt 


(a H f U J \V TV ^ v 

• [a"6“c“ff (?+ pt) + a-b-c+H (*+, p^ p~ z ) + (*>+ ^ , <flt ) 

+a~b^c ¥ H (px } py , + ^ + 6 c (<£ x ? ipt) + a ^ c 

-fa + ?> + c~i7 (<£~, + a + 6 + c + if (v? z :P y iP z )] 

b + b~ 


') (b + + b~) (c + + c~) 


(3.11) 


+ 


a a 


(o + + a - ) 


(pt - Px ) + 


(6+ + 6~) 


« - ^) + (c + + c -) ^ ' 


The three-dimensional scheme (3.11) can be further simplified by replacing a i j : ,a^j by 
di 3 and b ^ , b ~ 3 by bij similarly to the two-dimensional case, and also replacing j 

by aj = max± \{H Z pf)}\* In this case 

-I [#(</£ , (pU + H(pi, Py ’ ¥>7) + H(pt’Py ’ ¥?z ) + H(pt >py » ¥>z ) 
+#(<p~, <Pj) + H((pZ, pi,Pi) + H(Px’Py ’Pt) + ’Py-Pz )} ' 


(3.12) 
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3.3 A Dimension-by-Dimension 5th-order Reconstruction 

The reconstructions V p^ + p can be easily computed in a direction-by-direction manner. 
Such a direction- by-direction reconstruction is commonly used in upwind schemes [15], 
and we have used this strategy with central schemes in [6]. Here we show a three- 
dimensional example; generalizing this technique to more dimensions is straightforward. 
Using the notation of Section 2.2, a three-dimensional fifth-order reconstruction is 


• For each j, k: pf = reconstruct.^' (±, p*j,k) 

• For each i , k\ pf = reconstruct jp' (±, pi,*,k) 

• For each i,j: <pf = reconstruct.^' (±, 


where the subscript denotes the full range of an index: P*,j,k denotes the arraj 
(Vb j k, • • ■ > etc. We denote this operation in three dimensions as 




= reconstruct (±, p ) . 


(3.13) 




The results of this operation are three-dimensional arrays with elements ) iJk , 
)). k and (pf)ijj f 

Using this notation, we can turn Algorithm 2.1 into a three-dimensional scheme, 
replace reconstruct.^' (±, p n ) with reconstruct_V<p (±, p), and let F denotes the right 
hand side of (3.11). Applying this modified version of Algorithm 2.1 to each grid node 
gives a three-dimensional scheme that is fifth-order in space and fourth-order in time. 


4 Numerical Simulations 

In this section we present simulations that demonstrate the features of the schemes we 
developed in the previous sections. The scheme we test is the fifth-order semi-discrete 
method in one (Section 4.1), two (Section 4.2), and three (Section 4.4) space dimensions. 
Some of these examples are standard test cases that can be found, e.g., in [20, 28, 32]. 
In Section 4.3 we present a numerical stability study in two space dimensions. 

We do not follow the practice in [15] of masking singular regions from our error 
measurements, as we prefer to include the entire domain in our error estimate. 

4.1 One-Dimensional Examples 

A convex Hamiltonian 

We start by testing the performance of our schemes in a convex problem. We approxi- 
mate solutions of the one-dimensional equation 

d> t + ^ r + l) 2 = 0, ( 41 ) 

subject to the initial data <p(x ? 0) = — cos(7rx) with periodic boundary conditions on 
[0,2]. The change of variables, u(x,t) — c b x (x,t) + 1, transforms the equation into 
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the Burgers’ equation, u t + ± (u 2 ) x = 0, which can be easily solved via the method of 
characteristics [32], The solution develops a singularity in the form of a discontinuous 

derivative at time t = 7r . f 

The results of our simulations are shown in Figure 4.1. The order of accuracy o 
these methods is determined from the relative L 1 error, defined as the L -norm of t e 
error divided by the I^norm of the exact solution. These results along with the relative 
L~- norm before the singularity at T = 0.8/tt 2 , and after the singularity at T - 1.5/tt 
are given in Table 4.1. 


N 

— ^2 — 

Before singularity T = 0.8/7T 

relative L l -error 

L l -order 

relative L x -error 

L -order 

100 

2.78xl0 _b 

- 

5.74x10“' 

- 

200 

9.89xl0 _s 

4.81 

1.14x10“ 1 

5.65 

400 

3.20 xl0“ 9 

4.95 

1.92xl0“ lu 

5.90 

800 

1.01xl0 _i ° 

4.99 

3.04x10“^ 

5.98 

N 

After singularity T = 1.5/tt 2 

relative L l - error 

L l -order 

relative L°° -error 

L - order 

100 

2.04x 10 -4 

- 

2.02xl0“ 4 

— 

200 

7.21xl0 _? 

9.60 

8.15x10“' 

9.60 

400 

2.64 xlO" 3 

-5.19 

-5.19x10“° 

-6.65 

800 

2.55xl0 -5 

0.05 

0.05xl0“ b 

0.10 


Table 4.1: Relative L 1 -errors for the one-dimensional convex HJ problem (4.1) before 
(T = 0.8/tt 2 ) and after (T = 1.5/tt 2 ) the singularity formation. 


A non-convex Hamiltonian 

In this example we deal with non-convex Hamilton- Jacobi equations. In one dimension 
we solve 

(j) t — cos ( 4>x + 1) = 0, 

subject to the initial data <p (x, 0) = - cos (ttx) with periodic boundary conditions on 
[0 21 In this case (4.2) has a smooth solution for t < 1.049/tt , after which a singularity 
forms. A second singularity forms at t « 1-29/tt 2 . The results are shown in Figures 4.2. 
The convergence results before and after the singularity formation are given m lable 4.2. 

Remark. Tables 4.1 and 4.2 show that after the singularity formation the order of conver- 
gence deteriorates. In the following examples we will see that while a close examination 
of the convergence properties confirms this observation, in all the cases we examine , 
the error of the fifth-order semi-discrete scheme is less than the error of the other two 
published fifth-order methods for HJ equations [6, 15]. 
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Figure 4.1: One-dimensional convex Hamiltonian (4.1). Left: the solution before the 
singularity formation, T = 0.8/tt 2 . Right: the solution after the singularity formation, 
T = 1.5/7T 2 . N = 40. The fifth-order approximation is plotted on top of the exact 

solution. 


N 

Before singularity T = 0.8/n 2 

relative L l -error 

L 1 -order 

relative L°° -error 

L°°-order 

100 

1.20x10-* 

- 

4.24xl0 -7 

- 

200 

5. 29x10“* 

4.50 

2.18x10-* 

4.28 

400 

2.14xl0~ y 

4.62 

6.06xl0- 10 

5.17 

800 

8.24xl0 _il 

4.70 

1.17xl0 -11 

5.69 

N 

After singularity T = I.o/tt 2 

relative L l -error 

L l -order 

relative L°° -error 

L°°-order 

100 

1.91xl0' 5 

- 

3.52x10-* 

- 

200 

4.98xl0“ b 

-1.38 

4.27xl0- b 

-0.27 

400 

2.91xl0 _ti 

4.10 

2.47xl0- ti 

4.11 

800 

4.20xlO _B 

-0.53 

3.63x10-* 

-0.56 


Table 4.2: Relative L 1 -errors for the one-dimensional non-convex HJ problem (4.2) 
before (T = 0.8/tt 2 ) and after (T = 1.5/7T 2 ) the singularity formation. 
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Figure 4 2- One-dimensional non-convex Hamiltonian (4.2). Left: The solution before 
the singularity formation, T = 0.8/tt 2 . Right: The solution after the singularity forma- 
tion, T = 1.5/7T 2 . N = 40. The fifth-order approximation is plotted on top of the exact- 

solution. 

4.1.1 A comparison with existing fifth-order WENO-based methods 

In Figure 4.3 we compare the error of our new fifth-order semi-discrete scheme (Algo- 
rithm 2 1) with our fully-discrete scheme [6], and with the upwind WENO method of 
[15] (with a local Lax-Friedrichs flux). We also present results obtained with the method 
of [15] where the local Lax-Friedrichs flux was replaced by the semi-discrete central flux 
(2 7) which compares our smoothness measures with those of [15]. 

We see that before the singularity formation the terror of our semi-discrete method 
is as much as an order of magnitude smaller than the Zd-error of the methods m [6] and 
[151. The method of [15] with the flux (2.7) yields somewhat smaller errors for large grid 
spacing for the convex Hamiltonian, but becomes comparable to our method as the grid 
spacing decreases. For the non-convex Hamiltonian the method of [15] with flux (2.7) 
has larger errors than Algorithm 2.1. We take this as indication that the smoothness 
measures in [15] may be slightly better for large grid spacing and some Hamiltonians. 

After the formation of the singularity the behavior of the error in both methods that 
are based on the flux (2.7) is more erratic than the other two methods. Nonetheless, 
the methods that use the flux (2.7) have errors that are sometimes dramatically smaller 
than the other two methods. At no time is the error of methods using the flux (2.<) 
larger than that of the other two methods. Further comparisons are done in the next 
example and in Section 4.3. A theoretical study of the convergence of these schemes is 
beyond the scope of this work and is left for the future. 


A linear advection equation 

In this example ([15] with a misprint, corrected in [34]) we solve the one-dimensional 
linear advection equation, i.e., H {4> x ) = 4> x . We assume periodic boundary conditions 
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Figure 4.3: Convergence results for the convex Hamiltonian (4.1) (top) and non-convex 
Hamiltonian (4.2) (bottom). The relative iT-error are plotted against the number of grid 
nodes, “x": our semi-discrete fifth-order method (Algorithm 2.1). the fifth-order 

method of [6]. “o”: the fifth-order method of [15] with a local Lax-Friedrichs flux, “o”: 
the fifth-order method of [15] with the flux (2.7). Left: Before the singularity. Right: 
After the singularity. 
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on [—1, 1], and take the initial data as <fi (x, 0) = g (x — 0.5) on [ 1,1], where 
g (x) = - ^ {x + l) + h(x), 

( 2 cos (y x 2 ) — \/3, — l<x<— 5, 

,, x _ 1 3/2 + 3 cos(2to) , -|<x<0, (4.3) 

| 15/2 - 3cos (27 tt) , 0<x<5, 

[ (28 + 4tt + cos (37 tx)) /3 + 67 tx (x - 1) , 5 < x < 1. 

The results of our semi-discrete fifth-order method (Algorithm 2.1) are shown in Fig- 
ure 4.4, where it is compared with the fifth-order methods of [6] and [15]. The semi- 
discrete method (Algorithm 2.1) shows reduced dissipation compared to the method in 
[15], In [6] we showed that the fully-discrete fifth-order method we developed there is 
more stable than the method of [15] from the point of view of being able to use larger 
time steps. The numerical results here are based on fitting to each scheme its optimal 
time-step, hence the reduced dissipation for the fully-discrete scheme [6]. 


4.2 Two-Dimensional Examples 

A convex Hamiltonian 

In two dimensions we solve a problem similar to (4.1) 

<j>t + 2 (0* + fiy + l) 2 = 0, (4-4) 

which can be reduced to a one-dimensional problem via the coordinate transformation 
( f \ = ( V 2 1 / 2 1 ( X 1 The results of the second-order calculations for the 

l rj ) V 1/2 -1/2 ) \ V J 

initial data <f>(x,y,0) = — cos(7r(x + y)/ 2) = — cos(7r£) are shown in Figure 4.5. e 
convergence rates for the two-dimensional fifth-order scheme (3.9) before and after the 
singularity are shown in Table 4.3. 

A non-convex Hamiltonian 

The tw'o-dimensional non-convex problem, which is analogous to the one-dimensional 
problem (4.2), is 

<pt — cos ( 4>x + <py "F 1 ) — 0- (4-- 1 ) 

We assume the initial data <j> (x, y, 0) = - cos (tv(x + y)j 2), and periodic boundary con- 
ditions. The results are shown in Figure 4.6. The convergence results for the two- 
dimensional fifth-order scheme (3.9) before and after the singularity formation are given 
in Table 4.4. 
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■ — - — 2 ~ 

Before singularity T — 0.8/ 7r 

N 

relative L l -error 

L l -order 

relative L°° -error 

L°° -order 

50 

3.38xl0~ b 

- 

3.66xl0~ 7 

- 

100 

1.90xl0 _b 

4.15 

5.30xl0" g 

6.11 

200 

7.35x 10 -8 

4.69 

6.02xl0" u 

6.46 


After singularity T = 1.5/7T 

N 

relative L l -error 

L l -order 

relative L x -error 

L -order |; 

50 

8.68xl0“ 4 

- 

1.60xl0 _t> 

| 

100 

3.06xl0“ 4 

1.50 

2.88x10^ 

2.47 

200 

2.77xl0 -5 

3.46 

8.29xl0 _b 

5.12 "j 


Table 4.3: Relative L 1 - and L°°-errors for the two-dimensional convex HJ problem (4.4) 
before and after singularity formation, computed with (3.9) integrated in time as in 
Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 



Figure 4.5: Two-dimensional convex Hamiltonian, (4.4). Left: the solution before the 
singularity formation, T = 0.8/ 7r 2 . Right: the solution after the singularity formation, 
T = 1.5/7 r 2 . N = 40 x 40. The solution is computed with (3.9) integrated in time as in 
Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 
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Before singularity T = 0.8/ir 2 

N 

relative L l -error 

L l -order 

relative L 00 -error 

L°°-order 

50 

1.70xl0 _& 

- 

6.04xl0~ e 

- 

100 

1.69x 10 -ti 

3.33 

5.20x10-^ 

3.54 

200 

8.16xl0~ s 

4.37 

1.17xlO -iU 

5.47 


After singularity T = 1 .5/tt 

N 

relative L 1 -error 

L l -order 

relative L°°-error 

L°° -order 

50 

2.63xl0 -3 

- 

9.55xl0 _b 

- 

100 

3.40xl0~ 4 

2.95 

1.52x 10~ b 

2.65 

200 

7.20xl0 -& 

2.24 

2.72xl0- v 

2.49 


Table 4.4: Relative L 1 - and L°°-errors for the two-dimensional non-convex HJ problem 
(4.5) before and after the singularity formation, computed with (3.9) integrated in time 
as in Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 



Figure 4.6: Two-dimensional non-convex Hamiltonian, (4.4). Left: the solution before 
the singularity formation, T = 0.8/tt 2 . Right: the solution after the singularity forma- 
tion. T = 1 .5 / tt 2 . N = 40 x 40. The solution is computed with (3.9) integrated in time 
as in Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 
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A fully two-dimensional example 

To check the performance of our method on fully two-dimensional problems we solve a 
test problem which we introduced in [4]: 

<k + 4>x<i>y = 0, (x, y) € [-7T, 7T] X [-7T, 7r] , (4-6) 

subject to the initial data <f> (x, y, 0) = sin (x) + cos (y) and to periodic boundary 
conditions. The exact solution for this problem is given implicitly by 4>{x,y,t) = 
— cos ( q ) sin (r) + sin ( q ) + cos (r) where x = q — t sin (r) and y = r + t cos (q). This 
solution is smooth for t < 1, continuous for all t and has discontinuous derivatives for 
t > 1. The results of our simulations at times T = 0.8, 1.5, are shown in Figure 4.7. 
For comparison we show in Figure 4.8 the results obtained for the same problem with 
our fully-discrete method [6]. The convergence results for the fifth-order method (3.9) 
before the singularity formation are given in Table 4.5 and confirm the expected order 
of accuracy. 


Before singularity T = 0.8 

N 

relative L l -error 

L l -order 

relative L°° -error 

L°° -order 

50 

2.39xl0~ b 

- 

1.34x10-* 

- 

100 

8.52x10-* 

4.81 

1.40xl0 -11 * 

6.57 

200 

3.05x10-** 

4.80 

1.24X10" 12 

6.83 


Table 4.5: Relative L 1 -errors for the two-dimensional HJ problem (4.6) before the sin- 
gularity formation. T — 0.8. The solution is computed with (3.9) integrated in time as 
in Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 


An eikonal equation in geometric optics 

We consider a two-dimensional non-convex problem that arises in geometric optics [17] 

/ <f>t + + = 0’ (4.7) 

( <j> (x, y, 0) = j (cos (27 rx) — 1) (cos (27 ry) — 1) — 1. 

The results of our fifth-order method at time T = 0.6 are shown in Figure 4.9, where 
we see the sharp corners that develop in this problem. 

An optimal control problem 

We solve an optimal control problem related to cost determination [32]. Here the Hamil- 
tonian is of the form H(x, y, 

j (p t - sin (y) <t> x + sin (x) <f> y + \<f> y \ - \ sin 2 (y) - 1 + cos (x) = 0, r 4i g) 

1 4>(x,y, 0 ) = 0. 

The result of our fifth-order semi-discrete scheme at time T = 1 is shown in Figure 4.10. 
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Figure 4.7: Fully two-dimensional Hamiltonian, (4.6). Left: the solution before the 
singularity formation, T = 0.8. Right: the solution after the singularity formation, 
r = 1.5. N = 50 x 50. The solution is computed with (3.9) integrated in time as in 
Algorithm 2.1, with the fifth-order reconstruction of Section 3.3. 



Figure 4.8: The fully two-dimensional Hamiltonian computed with the method in [6] 
after the singularity formation, T = 1.5. N = 50 x 50. 
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Figure 4.9: Two-dimensional eikonal equation, (4.7). N = 40 x 40. Left: the initial 
data. Right: our semi-discrete fifth-order semi-discrete approximation at T = 0.6. 



Figure 4.10: Two-dimensional optimal control problem, (4.8). An approximation with 
our semi-discrete fifth-order method (3.9) is shown at T = 1. N = 40 x 40. 
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4.3 A Stability Study 

In this section we present a stability study, checking the stability properties of the 
two-dimensional semi-discrete fifth-order method. We compute the relative L l errors 
for various examples while varying the CFL number. In Figure 4.11 we compare the 
results obtained with our fifth-order scheme with the fully-discrete method [6], and 
with the upwind method of [15] using a local Lax- Friedrichs flux. As expected, the 
stability properties of the method (3.9) are similar to the stability properties of the 
upwind WENO method of [15], though our new method (3.9) enjoys smaller L 1 errors 
and hence is more accurate. 


4.4 Three-Dimensional Examples 

Finally, we solve a couple of three-dimensional problems with the scheme (3.11) inte- 
grated in time as in Algorithm 2.1, with the fifth-order reconstruction (3.13). We start 
with a convex problem 

<f>t + ~ (0* + + l) 2 = 0) 

subject to the initial data z, 0) = — cos (7r(x + y + z)/3). The convergence results 

for the scheme (3.11) before and after the singularity formation are given in Table 4.6. 
We also use (3.11) to approximate the solution of the non-covex problem 


0 t - cos (<f> x + <t> y + 4> z + 1) = 0, 


(4.10) 


with the same initial data. The convergence rates are shown in Table 4.7. 


Before singularity T = 0.5/7T 2 

N 

relative L l -error 

L l -order 

relative L°°-error 

L°°- order 

25 

1.04xl0~ 4 

- 

3.10xl0“ 8 

- 

50 

6.52xl0 _tl 

3.99 

2.66xl0 -lu 

6.87 

O 
o 
1 — 1 

3.74xl0" 7 

4.12 

2.02x10'^ 

7.04 

After singularity T — 1.5/7T 2 

N 

relative L l -error 

L 1 -order 

relative L 00 -error 

L°° -order 

25 

1.40xl0" 3 

- 

9.76xl0 -6 

- 

50 

1.80xl0 -4 

2.95 

4.15xlO~ B 

1.23 

100 

1.26xl0 -4 

0.51 

6.94x10-' 

2.58 


Table 4.6: Relative L l - and L°°-errors for the three-dimensional convex HJ problem 
(4.9) before and after the singularity formation, computed with (3.11) integrated in 
time as in Algorithm 2.1, with the fifth-order reconstruction (3.13). 
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Figure 4.11: Stability of the two-dimensional semi-discrete methods. A 1 — 100 x 100. 
“x”: our semi-discrete fifth-order method (3.9). our fully-discrete fifth-order 

method [6j. “o”: the fifth-order upwind method of [15] with a local Lax-Friedrichs 

flux. Upper left: linear advection with initial condition <j> (. x , y, 0) = - cos (n(x + y)/ 2). 
Upper right: fully 2D Hamiltonian (4.6). Middle row: convex Hamiltonian (4.4), be- 
fore the singularity (left) and after the singularity (right). Bottom row: non-convex 
Hamiltonian (4.5), before the singularity (left) and after the singularity (right). 
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^ » ■ n 

Before singularity T = 0.5/7T 

N 

relative L 1 -error 

L l -order 

relative L°°-error 

L°° -order 

25 

9.10xl0 _5 

- 

2.58x10-* 

- 

50 

3.85x10"“ 

4.56 

2.27xlO" IU 

6.83 

100 

1.77x 10~ 7 

4.45 

1.53xlO" rJ 

7.21 

After singularity T = 1.5/ 7T 2 

iV 

relative L l -error 

L 1 -order 

relative L°°-error 

L°° -order 

25 

9.99xl0 -4 

- 

6.60xl0" 7 

- 

50 

1.09xl0 -4 

3.20 

5.25xl0~ 7 

0.33 

100 

1.07x10"^ 

3.34 

6.13x10"“ 

3.01 


Table 4.7: Relative L l - and L°°-errors for the three-dimensional non-convex HJ problem 
(4.10) before and after the singularity formation, computed with (3.11) integrated in 
time as in Algorithm 2.1, with the fifth-order reconstruction (3.13). 


Appendix A: A Proof of Theorem 2.1 

Proof. Let H ( u ) 6 C 2 be convex (H"(u) > 0 or H" ( u ) < 0). We need to show that 
the flux 

H KKP (u\u-) = [a-H (Up + a+H („-)] - (u + - u) , 

is a non-increasing function of and a non-decreasing function of u~ . Here a + and 
a~ are defined as 


a + = max {H f (u) , 0} , a = min |{/T (u) , 0} | , 

u£l(u~ ,u Jr ) 


where / (a, b ) is the closed interval with endpoints a and b. 

The proof for u 4 * is discussed in detail. The proof for u~ is similar. 
Let u f > u- 2 ' Define the difference 


D = H KNP (ut,u-)-H KNP (i4,u~) 
1 


af + Oj 


H (uf) + af H (u )] 


a{ 

af + 


«- 


u 


1 


0-2 + a 2 


'a^H{ut)+a+H(u-)} + 




aj a 


2 u 2 


0,2 “ I - 0-2 


\Un 


We will prove that D < 0. We rewrite the D as the difference 


-) 

u ~ ) • 


D = G (ut) - G (u $) , 


where for fixed u 

G (u) = A (u) [H ( u ) - H (u~) - a + (u) (« - «")] , 
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with 


A(u) 


a ( u ) 

a + ( u ) + or (u) ’ 


a + (ti) 


max {H f (tz) ,0}, a (u) = 

u€l(u~ ,u ) 


min \{H* (u) , 0}| . 

,u) 


Since U+ > itj, the requirement D < 0 is equivalent to G' ( it ) < 0. Because H is 
convex (so the extrema of H' on an interval occur at the endpoints of that interval), 
a ± € C l so 


G'{u) = A' (it) [H (it) - H (tT) - a + (u) (it - «")] 

+A ( it ) [ H ' ( it ) — a +l (it) (u — u~) — a + (u)] 

Because H is continuous, there exists a £ (it) € I (u~,u) such that 
H'{£(u)){u-u-) = H(u)-H(u-) so 

G'{u) = A!(u){H'(i(u))-aA(u)] (u - u~) - A (u) a +l (u) (u - u~) 
+A ( it ) (H' ( it ) — a + (it)) 

= B ( it ) ( u — u ~ ) + A (tt) ( H ' ( u ) — a + (it)) , 

where B ( it ) = A' ( tt ) (H‘ (£ ( u )) — a + ( u )) — A (it) a +r (it) . 

Now 

a+ (it) or' (it) - a~ (it) a +/ (it) 

(o+ (it) + a" (u)) 2 

hence 

B(tt) = [G + (u)a-'(n)(^'(eH)-a + (^)) 

(a+ (u) + a - (it)) 

-a" (it) a + ' (u) (#' (£ (it)) + a" (it))] . 


We are now in a position to prove that D < 0. There are three cases to consider. 

Case 1. ut > tt 2 > it". In this case u - u~ > 0, and if ixi > ti 2 then [it~,iti] D [tt~,u 2 ] 
so o ± (iti) > a± (n 2 ) and aA' ( it ) > 0. Then 


B(u) 


a + (it 


)o-'(«) (H'(a«))-a + («)) 


> 0 


(a + ( u ) + a~ (it))' 
a~ (it) aA' (it) (H' (f (it)) + a - (it)) 


< 0 


> 0 


> 0 


< o, 


so 

B(tt)(tt — u _ ) A (u) (H' (u) - a + (u)) 

G' (it) = v ' + '<0. 

<0 >0 >0 <0 


Case 2. u~ > uf > itj. In this case it — it <0, and if iti > it 2 then [it , iti] C [it , u 2 ] 
so a ± (iti) < a ± (it 2 ) and a ±/ (it) < 0. We therefore have B (it) > 0 and 


B (it) (it -it ) H(u) (H' (it) - a + (it)) 

>0 <0 >0 <0 


G' (it) = 
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Case 3. uf >u~ >u%. In this case the proof is much more direct: by the continuity 
of H there exists a 6 € and a 6 [u^,u~] such that 


D = 


tti + CL] 


(H 1 (6) - at) K - u ~) ~ ( H ' (&) - 4 ) K - «•) 

v do ~r 0^2 


< o. 

The proof that H KNP (u + ,u~) is non-decreasing in u - is the same with 

G ( u ) = TTTT^ r ~\ \- H ( u ) " H ( u+ ) + a_ ^ ( U ~ 

v ' a + (u) + a~ (u) 

for fixed 


References 

[1] Abgrall R., Numerical discretization of the first-order Hamilton- Jacobi equation on 
triangular meshes , Comm. Pure Appl. Math., 49 (1996), pp.1339 1373. 

[2] Barles G., Solution de viscosite des equations de Hamilton- Jacobi, Springer- Verlag, 
Berlin, 1994. 

[3] Bianco F., Puppo G., Russo G., High order central schemes for hyperbolic systems 
of conservation laws, SIAM J. Sci. Comp., 21 (1999), pp. 294— 322. 

[4] Bryson S., Levy D., Central schemes for multi-dimensional Hamilton- Jacobi equa- 
tions, NAS Technical Report NAS-01-014, 2001, SIAM J. Sci. Comp., to appear. 

[5] Bryson S., Levy D., High-order central WE NO schemes for ID Hamilton- Jacobi 
equations, Proc. Enumath 2001, Ischia, Italy, to appear. 

[6] Bryson S.. Levy D., High-order central WENO schemes for multi- dimensional 
Hamilton- Jacobi equations, NAS Technical Report NAS-02-004, submitted to SIAM 
J. Nurner. Anal. 

[7] Crandall M.G., Evans L.C., Lions P.-L., Some properties of viscosity solutions of 
Hamilton- Jacobi equations, Trans. Amer. Math. Soc., 282 (1984), pp. 487-502. 

[8] Crandall M.G., Ishii H., Lions P.-L., User’s guide to viscosity solutions of second 
order partial differential equations, Bull. Amer. Math. Soc., 27 (1992), pp.1-67. 

[9] Crandall M.G., Lions P.-L., Viscosity solutions of Hamilton- Jacobi equations, Trans. 
Amer. Math. Soc., 277 (1983), pp.1-42. 

[10] Crandall M.G., Lions P.-L., Two approximations of solutions of Hamilton- Jacobi 
equations, Math. Comp., 43 (1984), pp.1-19. 


High-Order Semi-Discrete Schemes for HJ 


33 


[11] Friedrichs K.O., Lax P.D., Systems of conservation equations with a convex exten- 
sion , Proc. Nat. Acad. Sci . , 68 (1971), pp. 1686-1688. 

[12] Gottlieb S., Shu C.-W., Tadmor E., Strong stability-preserving high order time 
discretization methods, SIAM Review, 43 (2001), pp. 89-1 12. 

[13] Harten A., Engquist B., Osher S., Chakravarthy S., Uniformly high order accurate 
essentially non- oscillatory schemes III, J. Comp. Phys., 71 (1987), pp.231-303. 

[14] Kruzkov S.N., The Cauchy problem in the large for nonlinear equations and for 
certain quasilinear systems of the first order with several variables, Soviet Math. 
Dokl., 5 (1964), pp. 493-496. 

[15] Jiang G.-S., Peng D., Weighted ENO schemes for Hamilton- Jacobi equations, SIAM 
J. Sci. Comp., 21 (2000), pp.2126-2143. 

[16] Jiang G.-S., Shu C.-W., Efficient implementation of weighted ENO schemes, J. 
Comp. Phys., 126 (1996), pp.202-228. 

[17] Jin S., Xin Z., Numerical passage from systems of conservation laws to Hamilton- 
Jacobi equations and relaxation schemes, SIAM J. Numer. Anal., 35 (1998), pp.2385- 
2404. 

[18] Kurganov A., Lev>' D., A Third-Order Semi-Discrete Scheme for Conservation 
Laws and Convection- Diffusion Equations, SIAM J. Sci. Comp., 22, (2000), pp.1461- 
1488. 

[19] Kurganov A., Noelle S., Petrova G., Semi-discrete central-upwind schemes for hy- 
perbolic conservation laws and Hamilton- Jacobi equations, SIAM J. Sci. Comp., 23 
(2001), pp. 707-740. 

[20] Kurganov A., Tadmor E., New high-resolution semi-discrete central schemes for 
Hamilton- Jacobi equations, J. Comp. Phys., 160 (2000), pp.i20— <42. 

[21] Kurganov A., Tadmor E., New high-resolution central schemes for nonlinear conser- 
vation laws and convection- diffusion equations, J. Comp. Phys., 160 (2000), pp.241- 
282. 

[22] Levy D., Puppo G., Russo G., A fourth order central WENO scheme for multi- 
dimensional hyperbolic systems of conservation laws, SIAM J. Sci. Comp., to appear. 

[23] Levy D., Puppo G., Russo G., Central WENO schemes for hyperbolic systems of 
conservation laws, Math. Model, and Numer. Anal., 33, no. 3 (1999), pp.547 571. 

[24] Levy' D.. Puppo G., Russo G., Compact central WENO schemes for multidimen- 
sional conservation laws, SIAM J. Sci. Comp., 22 (2000), pp. 656-672. 

[25] Lions P.L., Generalized solutions of Hamilton- Jacobi equations, Pitman, London, 
1982. 



34 


S. Bryson and D. Levy 


[26] Lions P.L., Souganidis P.E., Convergence of MUSCL and filtered schemes for scalar 
conservation laws and Hamilton- Jacobi equations , Numer. Math., 69 (1995), PP-441- 
470. 

[27] Lin C.-T., Tadmor E., L 1 -stability and error estimates for approximate Hamilton- 
Jacobi solutions, Numer. Math., 87 (2001), pp. 701-735. 

[28] Lin C.-T., Tadmor E., High-resolution non- oscillatory central schemes for approx- 
imate Hamilton- Jacobi equations, SIAM J. Sci. Comp., 21, no. 6 (2000), pp.2163- 
2186. 

[29] Liu X.-D., Osher S., Chan T., Weighted essentially non- oscillatory schemes, J. 
Comp. Phvs., 115 (1994), pp. 200-212. 

[30] Nessyahu H., Tadmor E., Non- oscillatory central differencing for hyperbolic conser- 
vation laws, J. Comp. Phys., 87, no. 2 (1990), pp. 408-463. 

[31] Osher S., Sethian J., Fronts propagating with curvature dependent speed: algorithms 
based on Hamilton- Jacobi formulations, J. Comp. Phys., 79 (1988), pp. 12-49. 

[32] Osher S., Shu C.-W., High-order essentially nonoscillatory schemes for Hamilton- 
Jacobi equations, SIAM J. Numer. Anal., 28 (1991), pp.907-922. 

[33] Souganidis P.E., Approximation schemes for viscosity solutions of Hamilton- Jacobi 
equations, J. Diff. Equations., 59 (1985), pp.1-43. 

[341 Zhang Y.-T., Shu C.-W., High-order WENO schemes for Hamilton- Jacobi equations 
on triangular meshes, NASA/CR-2001-211256, ICASE Report No. 2001-39, (2001), 
submitted to SLAM J. Sci. Comp. 



